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We put forward a simpler and improved variation of a recently proposed method to overcome the 
signal-to-noise problem found in Monte Carlo calculations of the entanglement entropy of interacting 
fermions. The present method takes advantage of the approximate lognormal distributions that 
characterize the signal-to-noise properties of other approaches. In addition, we show that a simple 
rewriting of the formalism allows circumvention of the inversion of the restricted one-body density 
matrix in the calculation of the n-th Renyi entanglement entropy for n > 2. We test our technique 
by implementing it in combination with the hybrid Monte Carlo algorithm and calculating the 
n = 2, 3,4,..., 10 Renyi entropies of the ID attractive Hubbard model. We use that data to 
extrapolate to the von Neumann (n = 1) and n —>■ oo cases. 
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INTRODUCTION 

Recently [1], we proposed an algorithm to compute the 
Renyi entanglement entropy Sn of interacting fermions. 
Many algorithms have been proposed to this effect in 
the last few years [2-8]. Our proposal, based on the 
free-fermion decomposition approach of Ref. [10], over¬ 
comes the signal-to-noise problem present in that ap¬ 
proach and is compatible with the hybrid Monte Carlo 
(HMC) method [9] widely used in the context of lattice 
quantum chromodynamics. The core idea of our method 
is that, by differentiating with respect to an auxiliary 
parameter A, one may carry out a Monte Carlo (MC) 
calculation of dSn/dX with a probability measure that 
includes entanglement properties explicitly. [This was 
not the case in the approach of Ref. [10], where the prob¬ 
ability measure factored across auxiliary field replicas; we 
identified this as the cause of the signal-to-noise problem 
(see below)]. Once the MC calculation is done, integra¬ 
tion with respect to A returns the desired entanglement 
entropy relative to that of a noninteracting system (which 
is easily computed separately). 

In this work, we describe and implement a variation 
on that Monte Carlo algorithm which, while sharing the 
properties and core idea mentioned above, differs from 
it in two important ways; the new method, in fact, is 
different enough that we advocate its use over our origi¬ 
nal proposal. First, the new method takes advantage of 
the approximate lognormal shape of the underlying sta¬ 
tistical distributions of the fermion determinants, which 
we already noted in Ref. [1] and which we explain in de¬ 
tail below. Second, and more importantly, the present 
method is simpler than our original proposal: whereas in 
the latter the parameter A multiplied the coupling con¬ 
stant g (thus generating a rather involved set of terms 
upon differentiation of the fermion determinant), here A 
is coupled to the number of fermion species N^. As we 
show below, this choice not only simplifies the implemen¬ 
tation, but also exposes the central role of the logarithm 


of the fermion determinant in our calculation of Sn, and 
thus brings to bear the approximate lognormality prop¬ 
erty mentioned above. 

Below, we present the basic formalism, review the ev¬ 
idence for approximate lognormal distributions, and ex¬ 
plain our method. Besides the points mentioned above, 
in our calculations we have found the present method 
to be more numerically stable than its predecessor. We 
explain this in detail in our Results section. 

In addition to the new method, we show that it is 
possible to rewrite part of the formalism in order to by¬ 
pass the calculation of inverses of the restricted density 
matrix (see e.g. [I, 6, 7]) in the determination of Renyi 
entropies of order n > 2. To test our method, we com¬ 
puted the n = 2 Renyi entropy of the ID attractive Hub¬ 
bard model using the previous as well as the new for¬ 
malism, and checked that we obtained identical results. 
Going beyond the n = 2 case, we present results for the 
n = 2,3,4,..., 10 Renyi entropies and find that higher 
Renyi entropies display lower statistical uncertainty in 
MC calculations. 

BASIC FORMALISM 

As in our previous work, we set the stage by briefly 
presenting the formalism of Ref. [10]. The n-th Renyi 
entropy S'„ of a sub-system A of a given system is 

•S'™ = :j-^lntr(p]^), (1) 

where is the reduced density matrix of sub-system A. 
For a system with density matrix p, the reduced density 
matrix is defined via a partial trace over the Hilbert space 
corresponding to the complement A of our sub-system: 

Pa = (2) 

An auxiliary-field path-integral form for p^, from 
which can be computed using MC methods for a wide 
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variety of systems, was presented in Ref. [10], which we 
briefly review next. 

As is well known from conventional many-body formal¬ 
ism, the full density matrix p can be written as a path 
integral by means of a Hubbard-Stratonovich auxiliary- 
field transformation: 

g-/3A r 

P = = / 'DaPla] p[ct], (3) 


an explicit form can be identified in the numerator as in 
the replica trick [11], which corresponds to a partition 
function for n copies of the system, “glued” together in 
the region A. 

The naive probability measure, namely 

n 2N 

p{w }]=n n (10) 

fc=l m=l 


for some normalized probability measure T’[cr] deter¬ 
mined by the details of the underlying Hamiltonian (for 
more detail, see below and also Ref. [12]). Here, Z is 
the partition function, and p[a\ is the density matrix of 
noninteracting particles in the external auxiliary held cr. 
One of the main contributions of Ref. [10] was to show 
that the above decomposition determines not only the 
full density matrix but also the restricted one. Indeed, 
Ref. [10] shows that 

Pa = j V(tP[<j]paW], ( 4 ) 

where P[a] is the same probability used in Eq. (3), 

PaW] = CaM exp l-J2cl[\n{GA^[a] - l)],^c^- 

\ ij 

(5) 

and 

Ca[c^] = det(l - Ga[<j])- (6) 

Here, G^[ct] is the restricted Green’s function of the 
noninteracting system in the external held cr (see below), 
and fil, c are the fermion creation and annihilation oper¬ 
ators. The sums in the exponent of Eq. (5) go over those 
points in the system that belong to the subsystem A. 

Using the above formalism for the case of 2N- 
component fermions, the entanglement entropy (c.f. 

Eq. 1) takes the form 

exp ((1 - n)S„) = J V{a}P[{a}] (3[{cr}], (7) 

where the held integration measure, given by 

= n (8) 

k=l 

is over the n “replicas” cr^ of the Hubbard-Stratonovich 
held (which result from taking the n-th power of the path 
integral representation of Pa shown above), and the nor¬ 
malization 

. 2N 

Z= VaUdetU^la] (9) 

m=l 

was included in the measure. It is worth noting that, by 
separating a factor of Z" in the denominator of Eq. (7), 



factorizes across replicas, which makes it insensitive to 
entanglement. This factorization is the main reason why 
using P[{cr}] as a MC probability leads to signal-to-noise 
issues (see Ref. [10]). In Eq. (10), [/^[ct] encodes the 
dynamics of the m-th fermion component, including the 
kinetic energy and the form of the interaction after a 
Hubbard-Stratonovich transformation. That matrix also 
encodes the form of the trial state |d>) in ground-state 
approaches (see e.g. Ref. [12]), which we use here; we 
have taken |d>) to be a Slater determinant. In finite- 
temperature approaches, Um[o'] is obtained by evolving 
a complete set of single-particle states in imaginary time. 

The quantity that contains the pivotal contributions 
to entanglement is 


2N 

QiM] = n detM„[{cr}], (11) 

m—1 

which we refer to below as the “entanglement determi¬ 
nant,” and where 


M„[{cr}] = (1 - G^_„[CTfc]) X 








( 12 ) 


The product Q[{cr}] played the role of an observable in 
Ref. [10], which is a natural interpretation given Eq. (7). 
However, we will interpret this differently below. Other 
than the field replicas, the new ingredient in the determi¬ 
nation of is the restricted Green’s function Ga ml'^k]- 
This is the same as the noninteracting one-body den¬ 
sity matrix G{x, x') of the m-th fermion component in 
the background field ct^., but the arguments x,x' are re¬ 
stricted to the region A (see Ref. [10] and also Ref. [13], 
where expressions were originally derived for the reduced 
density matrix of noninteracting systems, based on re¬ 
duced Green’s functions). 


AVOIDING INVERSION OF THE REDUCED 
GREEN’S FUNCTION FOR n > 2 

As noted in Ref. [14], for n = 2, no inversion of 
1 — Ga m i® actually required in the calculation of 
the entanglement determinant Q[{(t}], as the equations 
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clearly simplify in that case. However, for higher n it 
is not obvious how to avoid such an inversion. Here, 
however, we show that this calculation can indeed be ac¬ 
complished without inversion. We begin by noting that 


det M„[{cr}] = deX Lj,{a}]dei 


(13) 


where Ljn[{(^}] is a block diagonal matrix (one block per 
replica k): 


Lm[{cr}] = diag [1 - , 


(14) 


and 




with 


/ 1 

R[ai] 

0 

V 0 


0 0 ... 

1 0 ... 

i?[cr2] 1 0 


0 \ 

: 0 

0 0 


0 1 / 


R[<Jk] = 


G, 




(15) 


(16) 


The equivalence of the determinants in Eq. (13) can be 
shown in a straightforward fashion: the Lm[{(T}] factor is 
easily understood, as that matrix is block diagonal and 
therefore its determinant reproduces the first r.h.s. factor 
in the first line of Eq. (12); the remaining factor relies on 
the identity 


/a 0 0 ... 0 Hk\ 


det 


-Hi 1 0 ... 

0 -i?2 1 0 


= deX{t + HiH2...Hk) 


Vo .0 -Hk-i a ) 


(17) 

which is valid for arbitrary block matrices Hj, is a stan¬ 
dard result often used in many-body physics (especially 
when implementing a Hubbard-Stratonovich transforma¬ 
tion) , and can be shown using so-called elementary oper¬ 
ations on rows and columns. 

Within the determinant of Eq. (13), we may of course 
multiply Kjn[o'] and Lmicr]: 

Tm[{a}] = Km[{a}]L^[{a}] = a - ^^[{a}]^, (18) 

where Gm[{o'}] is a block diagonal matrix defined by 


SmiW}] = diag , 


(19) 


and 


B = 


/ a 0 
a a 
0 a 


... -a\ 
... 0 
... 0 


( 20 ) 


V 0 ... 0 a a / 


Equation (18) shows our claim, as we may use Tm[{<jY\ 
in our calculations instead of Mm[{<^W, and the former 
contains no inverses of a — G^ 

Summarizing, a class of approaches to calculating S'„ 
for n > 2, based on the Hubbard-Stratonovich represen¬ 
tation of Pa (also known as free-fermion decomposition), 
requires computing Mm[{(T}], which in turn requires in¬ 
verting a — Gji^ per Eq. (12). By arriving at Eq. (18), 
and given that 


detrm[{CT}] = detM„[{(T}], (21) 

[Eq. (13) and beyond] we have shown that no inversions 
are actually required, as Tmljc}] contains no inverses. 
While this is a desirable feature from a numerical point of 
view, it should be mentioned that, from a computational- 
cost point of view, the price of not inverting a — G^ ^ 
reappears in the fact that Tm, though sparse, scales lin¬ 
early with n in size. 

For the remaining of this work, calculations carried 
out at n = 2 use the M approach, which is based on 
Eq. (12) and the ‘proposed method’ described below. We 
reproduced those results by switching to the T approach, 
which uses Eq. (18) (as well as the method described 
below), and then proceeded to higher n with the latter. 


A STATISTICAL OBSERVATION: LOGNORMAL 
DISTRIBUTION OF THE ENTANGLEMENT 
DETERMINANT 


In Ref. [1], we presented examples of the approximate 
log-normal distributions obeyed by (5 [{(t}] when sampled 
according to H[{ct}]. One such example is reproduced 
here for reference in Fig. 1. The fact that such distri¬ 
butions are approximately log-normal, at least visually, 
suggests that one may use the cumulant expansion to 
determine S^. Indeed, in general. 


(I - n)S'„ = Iny V{a}P[{a}\ Q[{cr}] 

^ Pn Q] 

^ ml ’ 

m—1 


( 22 ) 


where K^pn Q] is the m-th cumulant of In Q, and the first 
two nonzero cumulants are given by 


«i[X] = (A) (23) 


and 


= (a2) - {Xf (24) 

for a functional A[{(7}], and where the expectation value 
(•) here taken with respect to the produce measure 
P[{(j}]. If the distribution of InQ were truly gaussian, 
the above series would terminate after the first two terms. 
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FIG. 1: (color online) Distribution of the observable Q)!®"}] 
of the naive free-fermion decomposition method, i.e. using 
Eq. (7), for a ten-site Hubbard model described by Eq. (33), 
at attractive coupling U/t = 2.0 and for a subsystem of size 
Lj^/L = 0.8. Here, Q[{cr}] is a non-negative quantity. The 
long tail in the main plot (note logarithmic scale in vertical 
axis) is approximately a log-normal distribution [i.e. In (3 [{(t}] 
is roughly a normal distribution (see inset)]. 


which would provide us with an efficient way to bypass 
signal-to-noise issues in the determination of with 
stochastic methods [15]. Unfortunately, the distribution 
is not exactly gaussian. Moreover, the cumulants be¬ 
yond m = 2 are often extremely sensitive to the details 
of the distribution (i.e. they can fluctuate wildly), they 
are hard to determine stochastically (the signal-to-noise 
problem re-emerges), and there is no easy way (that we 
know of) to obtain analytic insight into the large-m be¬ 
havior of K^. However, this approximate log-normality 
does provide a path forward, as it indicates that we may 
still evaluate (InQ) with good precision with MC meth¬ 
ods. As we will see in the next sections, this is enough 
to determine if we are willing to pay the price of a 
one-dimensional integration on a compact domain. 

Although (approximate) lognormality in the entangle¬ 
ment determinant seems very difficult to prove analyti¬ 
cally in the present case, evidence of its appearance has 
been found in systems as different as ultracold atoms and 
relativistic gauge theories [15, 16]. The underlying reason 
for this distribution appears to be connected to a similar¬ 
ity between the motion of electrons in disordered media 
and lattice fermions in the external auxiliary (gauge) held 
in MC calculations. 


PROPOSED METHOD 

Starting from the right-hand side of Eq. (7), we in¬ 
troduce an auxiliary parameter 0 < A < 1 and dehne a 


function r(A;( 7 ) via 


r(A;g)^ Jv{a}P[{a}] g^[W]. 

(25) 

At A = 0, 


lnr(0;5)=0. 

(26) 

while for A = 1, r(A; (/) yields the entanglement entropy: 

^\nT{yg) = S^. 

1 — n 

(27) 

Using Eq. (25), 


^=y V{a}P[{a}-X] lnQ[{a}] 

(28) 

where 


^[W;A] = ^^p[W] g^[W]. 

(29) 


In the presence of an even number of havors 2N and 
attractive interactions, T’[{cr}] and Q[{o'}] are real and 
non-negative for all a, such that there is no sign problem 
and P[{cr}; A] above is a well-dehned, normalized proba¬ 
bility measure. 

As in our previously proposed method, we can then 
calculate 5'„ by taking the A = 0 point as a reference and 
computing using 

S^ = -^ r dX{\nQ[{a}])^, (30) 

Jo 

where 

Wa = Jv{a}P{{ay,\]X[{a}]. (31) 

We thus obtain an integral form of the interacting Renyi 
entropy that can be computed using any MC method (see 
e.g. [12]), in particular HMC [9]. 

As in our previous work, we note that the above expec¬ 
tation values are determined with respect to the probabil¬ 
ity measure P[{(t};A] , which communicates correlations 
responsible for entanglement. In contrast to the canoni¬ 
cal MC probability P[{(t}], which corresponds to statis¬ 
tically independent copies of the Hubbard-Stratonovich 
held, this admittedly more complicated distribution does 
not exhibit the factorization to blame for the signal-to- 
noise problems present in the approach as originally for¬ 
mulated. 

Using Eq. (30) requires Monte Carlo methods to eval¬ 
uate (ln(5[{tT}])_,^ as a function of A, followed by integra¬ 
tion over A. As in our previous method, we hnd here that 
(In (5[{o'}])_)^ is a smooth function of A, which is essentially 
linear in the present case. It is therefore sufficient to per¬ 
form the numerical integration using a uniform grid. The 
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U/t = 0.5 
U/t = 1.0 
U/t = 2.0 


FIG. 2: (color online) Stochastic results for (In(5[{cr}])_n with 
n = 2 for couplings Ujt = 0.5,1.0, and 2.0 as functions of 
auxiliary parameter A and region size Lj^/L for a ten-site 
Hubbard model. 


stochastic evaluation of (In (^[{cr}])^,^, for fixed subregion 
A, can be expected to feature roughly symmetric fluctu¬ 
ations about the mean. As a consequence, the statistical 
effects on the entropy are reduced after integrating over 

A. 

Finally, we note an interesting application of Jensen’s 
inequality at A = 0. At that point 


d\nT 

dX 


A=0 


I V{a}P[{a}] lnQ[{a}] 


(32) 


< Iny V{a}P[{a}] Q[{a}] = (1 - n)S'„, 


which must be satisfied by our calculations. Our Monte 
Carlo results at A = 0 indeed satisfy this bound. 


RESULTS 

Second Renyi entropy 

As a first test of our algorithm and in efforts to make 
contact with previous work [1, 10], we begin by show¬ 
ing results for the second Renyi entropy S 2 for the one¬ 
dimensional Hubbard chain with periodic boundary con¬ 
ditions at half filling, whose Hamiltonian is 

^ = + CysCi.s) 

where the first sum includes s =t,4. and pairs of ad¬ 
jacent sites. We implemented a symmetric Trotter- 
Suzuki decomposition of the Boltzmann weight, with 
an imaginary-time discretization of r = 0.05 (in lattice 
units). As mentioned earlier, the many-body factor in 
the Trotter-Suzuki approximation was treated by intro¬ 
ducing a replica auxiliary field a for each power of the 


reduced density matrix. As in our previous work, we im¬ 
plemented a Hubbard-Stratonovich transformation of a 
compact continuous form [12]. 

We present plots for (lnQ[{o'}])A with n = 2 in Fig. 2. 
In contrast to the results obtained in Ref. [1] and as men¬ 
tioned above, the resulting expectation demonstrates sur¬ 
prisingly little curvature as the region size is varied 
and is stunningly linear as a function of the auxiliary 
parameter A. Even after twice doubling the strength of 
the interaction, the curvature of constant-subsystem-size 
slices is increased only marginally. We note that if one 
assumes such benign curvature is a somewhat universal 
feature, at least for weakly-coupled systems, our method 
provides a means by which to rapidly estimate the entan¬ 
glement entropy for a large portion of parameter space at 
the very least yielding a qualitative picture of its behavior 
as a function of the physically relevant input parameters. 

We observe that this surface displays almost no torsion, 
its dominant features being those present in the non¬ 
interacting case i.e. an alternating shell-like structure. 
Toward larger region sizes, we observe a combination of 
twisting and translation culminating in the required, and 
somewhat delicate, cancellation upon reaching the full 
system size. Presented with this relatively forgiving ge¬ 
ometry, we performed the required integration via cubic- 
spline interpolation. Using a uniformly spaced lattice of 
size = 20 points, we determine the desired entropy to 
a precision limited by statistical rather than systematic 
considerations. 



FIG. 3: (color online) Resnlts for the ten-site Hubbard chain 
for couplings U/t = 0.5,1.0, 2.0, and 4.0 for 7,500 samples 
with associated numerical uncertainties. Results for U/t = 0 
are included as a dashed line (black). For all but the largest 
coupling, exact diagonalization results from Ref. [10] are indi¬ 
cated by solid lines, while for the largest coupling, we provide 
a line joining the central values of our result to emphasize 
that its shape is consistent with results for the former. 
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FIG. 4: (color online) Entanglement entropy S 2 in nnits of the 
result for a free system plotted as a function of the number of 
samples Ns for couplings U/t = 0.5,1.0, 2.0, and 4.0 demon¬ 
strating convergence to within a few percent within the first 
ten thousand samples. 


Comparison to exact diagonalization 

Shown in Fig. 3 are results for a system of size L = 
with a number of sites = 10. For couplings U/t = 
0.5,1.0, 2.0 and 4.0 and region sizes = 1,2, ...,10, 
we find solid agreement with previous calculations in 
Refs. [1, 10], and as in the former, we observe convergence 
rather quickly with only O(IO^) decorrelated samples as 
can be seen in Fig. 4. Further, for large sample sizes 
Ns, we observe that the standard error in the entropy 
AS' 2 , computed from the envelope defined by the MC 
uncertainty in the source (ln(5[{(T}])_,^ for at each value 
in (L^, A)-space, scales asymptotically as AS '2 ~ 
up to minute corrections. 


Results for n yf 2 


In this section, we extend the results presented above 
to n = 3,4, 5,..., 10. In order to highlight the differences 


between n = 2 and n > 2, we show in Fig. 5 the Renyi 
entropies S„ for n = 2,3,4 (top to bottom) of the ID 
attractive Hubbard model, as obtained with our method 
and the reformulation of the fermion determinant shown 
in Eq. (18). 

As evident from the figure, increasing n leads to lower 
values of S„ at fixed subsystem size LaIL consistent 
with knowledge that the Renyi entropy is a nonincreasing 
function of its order. However, increasing n also ampli¬ 
fies the fluctuations as a function of La/L- Interestingly, 
the approach of our system to the large-n regime is quite 
rapid, and after only the first few orders, the difference 
between consecutive entropies is only marginal, most ob¬ 
viously so at weak coupling. We also observe that, as 
n is increased, the statistical fluctuations that define the 
error bars appear to be progressively more suppressed, 
which is particularly evident for the strongest coupling 
we studied, namely U/t = 4.0. 

At the level of the auxiliary function (lnQ[{(T}])_)^, we 
again see very predictable changes in the geometry of this 
surface as a function both arguments as shown in Fig. 6. 
With fixed coupling and particle content, increasing the 
Renyi order results in a tilting effect reminiscent of that 
seen previously with increasing coupling, but rather than 
being localized away from vanishing subsystem size, the 
change is much more global, affecting all subsystems in 
a qualitatively similar fashion and leaving each surface’s 
characteristic quasi-linearity in A intact. Although the 
shell-like structure present in this function’s depen¬ 
dence is amplified, this increased fluctuation affects the 
quality of the results negligibly at most, as again, the 
geometry remains amenable to fairly naive quadratures. 

With the data presented above, we would be remiss 
if we did not attempt an extrapolation not only to the 
limit of infinite Renyi order Soo, but also to the von Neu¬ 
mann entropy, despite knowledge of the formidable chal¬ 
lenges presented by these extrapolations, particularly in 
the case of the latter. The former limit provides a lower 
bound on all finite-order entropies, whereas the latter is 
of interest to a variety of disciplines and has proven diffi¬ 
cult to study. At fixed coupling and with the knowledge 
that the Renyi entropy is nonincreasing in the order, we 
found that our results at each fixed region size and at 
every studied coupling were well-characterized by expo¬ 
nential decays. 

Interestingly, the relative speed of this decay oscillates 
as a function of the region size as can be seen in Fig. 7. 
Regions corresponding to an even number of lattice sites 
demonstrate a much more sudden initial decay than do 
those regions comprised of an odd number of sites. This 
peculiar oscillation results in an inverted shell structure 
for the extrapolation to n = I, in contrast to the case 
where n —> 00 in which this feature is preserved. A rep¬ 
resentative example of this procedure is shown in Fig. 8. 







7 



1.8 



0.6 

0 0.2 0.4 0.6 0.8 1 

La/L 



FIG. 6: (color online) Stochastic results for (In Q[{f’'}])A with 
n = 2, 4, 6, and 8 (top to bottom) for a coupling of U/t = 2.0 
as functions of auxiliary parameter A and region size Lj^jL. 



FIG. 7: (color online) Renyi entropies Sn for n = 2,4,6, 
and 8 (top to bottom with error bars and colors matching 
those in Fig. 6) of the ID attractive Hubbard model, as a 
function of the subsystem size La/L. The solid black line 
shows extrapolation to n = 1. The dashed black line shows 
extrapolation to n —>■ (xj. Again, results are shown for U/t = 
2 . 0 . 


FIG. 5: (color online) Renyi entropies S„ for n = 2, 3,4 (top 
to bottom) of the ID attractive Hubbard model, as a function 
of the subsystem size La/L. In each plot, results are shown 
for several values of the attractive coupling U/t. 


SUMMARY AND CONCLUSIONS 

We have presented a method to compute the entan¬ 
glement entropy of interacting fermions which takes ad¬ 
vantage of an approximate log-normality property of 
the distribution of fermion determinants. The result¬ 


ing approach overcomes the signal-to-noise problem of 
naive methods, and is very close in its core idea to an¬ 
other method we proposed recently; both methods in¬ 
volve defining an auxiliary parameter A, differentiating, 
and then integrating to recover Sn after a MC calcula¬ 
tion. The order of the steps is important, as the dif¬ 
ferentiation with respect to A induces the appearance of 
entanglement-sensitive contributions in the MC proba¬ 
bility measure. Beyond those similarities, the present 
method has the distinct advantages of being simultane¬ 
ously simpler to formulate (algebraically as well as com¬ 
putationally) and of explicitly using the approximate log- 
normality property. Moreover, we have found that the A 












FIG. 8: (color online) Interpolation of the Renyi entropies Sn 
for n = 2,3,4,10 for a coupling oi Ujt = 2.0 given as 
functions of the auxiliary parameter A as well as the region 
size Lj^/L. An extrapolation to n = 1 (the von Neumann 
entropy) as well as to n —> oo are shown in solid and dashed 
lines respectively. 


integration step displays clearly more stable numerical 
behavior in the present approach than in its predeces¬ 
sor: it is approximately linear in the present case and 
markedly not so in the original incarnation. We there¬ 
fore strongly advocate using the present algorithm over 
the former. 

In addition to presenting an improved method, we have 
put forward a straightforward algebraic reformulation of 
the equations which, while exactly equivalent to the origi¬ 
nal formalism, avoids the numerical burden of computing 
inverses of restricted Green’s functions in the calculation 
of n-th order Renyi entropies for n > 2. This issue had 
been pointed out by us and others (see e.g. Ref. [14]) 
as an inconvenience, as it is perfectly possible for those 
matrices to be singular. 

As a test of our algorithm, we have presented results 
for the Renyi entropy S'„ of the half-filled ID Hubbard 
model with periodic boundary conditions. The present 
and old formalisms were used for calculations at n = 2, 
which matched exactly. The rewritten form based on 
Eq. (18) was then used to extend our computations to 
n = 3,4,..., 10, allowing us to attempt extrapolations in 
the Renyi order in both directions. 

Our results show that, with increasing Renyi order n, 
the value of Sn decreases for all La/L, and the fluctu¬ 
ations as a function of LaIL become more pronounced. 
Remarkably, the statistical MC fluctuations decrease as n 
is increased. Since the problem we set out to solve was in 
fact statistical in nature, our observations indicate that 


calculations for large systems and in higher dimensions 
will benefit from pursuing orders n > 2. 
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